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We present stable bright solitons built of coupled unstaggered and staggered components in a 
symmetric system of two discrete nonlinear Schrodinger (DNLS) equations with the attractive self- 
phase-modulation (SPM) nonlinearity, coupled by the repulsive cross-phase-modulation (XPM) in- 
teraction. These mixed modes are of a "symbiotic" type, as each component in isolation may only 
carry ordinary unstaggered solitons. The results are obtained in an analytical form, using the varia- 
tional and Thomas- Fermi approximations (VA and TFA), and the generalized Vakhitov-Kolokolov 
(VK) criterion for the evaluation of the stability. The analytical predictions are verified against 
numerical results. Almost all the symbiotic solitons are predicted by the VA quite accurately, and 
are stable. Close to a boundary of the existence region of the solitons (which may feature several 
connected branches), there are broad solitons which are not well approximated by the VA, and are 
unstable. 

PACS numbers: 42.65.Tg; 05.45. Yv; 63.20.Ry; 03.75.Lm 



I. INTRODUCTION 

Discrete nonlinear Schrodinger (DNLS) equations constitute a class of lattice models which comprise di- 
verse physical settings A straightforward realization of the DNLS equation in arrays of evanescently 
coupled optical waveguides was first proposed in Ref. @, and later demonstrated experimentally in a set 
of parallel semiconductor waveguides Multi-core nonlinear waveguiding systems have also been created 
in the form of optically- written virtual lattices in photorefractive materials [|[ , and as permanent structures 
written by laser pulses in bulk silica [f|. A thorough review of the nonlinear discrete optics, developed 
experimentally and theoretically in these and allied media, was given in Ref. @ . The DNLS equations find 
another important application in modeling the mean-field dynamics of Bose-Einstein condensates (BECs) 
loaded into deep optical-lattice potentials. In this case, it was demonstrated experimentally Q and theoret- 
ically Q that the periodic potential effectively splits the condensate into a set of droplets trapped in local 
potential wells, which are linearly coupled by tunneling of atoms across the separating potential barriers, 
DNLS equations being natural models for such quasi-discrete systems. 

Two fundamental types of discrete solitons supported by the DNLS equations with the self-repulsive and 
self-attractive on-site nonlinearity are localized modes of staggered and unstaggered types, respectively, i.e., 
ones with opposite signs of the lattice field at adjacent sites, or without the sign alternation pj. In the 
continuum limit, the unstaggered solitons correspond to regular ones, residing in the semi-infinite gap of the 
continual NLS equation, while the staggered solitons may be considered as counterparts of gap solitons, which 
exist in finite bandgaps of the spectrum induced by a periodic potential, in the case of the self-dcfocusing 
nonlinearity [Til ]. 

A natural generalization, which also finds many applications to optics and BEC, is represented by systems 
of coupled DNLS equations. In optics, the system models the co-propagation of two waves carried by 
different polarizations or wavelengths in the same waveguiding array, while in BEC the coupled equations 



describe a mixture of two condensates, which may represent either different hyperfine states of the same" 
atomic species, or two different kinds of atoms |l0| . Normally, two-component discrete solitons in such 
systems feature one type of the intrinsic structure, unstaggcrcd or staggered, in both components, because 
the signs of the self-phase- modulation (SPM) nonlinearity acting on each component, and of the cross-phasc- 
modulation (XPM) nonlinearity which couples them, are the same The objective of the present work 
is to introduce two-component discrete solitons of the mixed type, built as combinations of unstaggered 
and staggered components. Previously, single-component surface modes of a mixed unstaggered-staggered 
type were studied at an interface between two different lattices [l2j|, but, to the best of our knowledge, 
two-component mixed solitons were not reported before. On the other hand, in bimodal continual systems 
with the periodic potential acting on both components, solitons of a semi-gap type, which may be considered 
as continuous counterparts of the discrete ones introduced in the present work, were studied in Ref. [l3|. 
They are composed of an ordinary soliton in one component and a gap soliton in the other. The semi-gap 
solitons are somewhat similar to the earlier studied intergap solitons, that were built as bound states of two 
components represented by solitons belonging to two different finite bandgaps (the first and second ones) 

The paper is organized as follows. The model is formulated in Section II. Approximate analytical results 
are presented in Section III. These results, based chiefly on the variational approximation (VA), demonstrate 
that the mixed unstaggered-staggered solitons are possible in the symmetric system of DNLS equations when 
the XPM interaction between the two components is repulsive, on the contrary to the self-attractive SPM 
nonlinearity. The situation with the opposite signs of the SPM and XPM terms seems exotic in optics, but it 
is quite possible in BEC, where the sign of the interactions may be readily switched by means of the Fcshbach 
resonance (see, e.g., Ref. [HI). In the case of strong difference between masses of the two components, another 
analytical solution is obtained, based on the Thomas- Fermi approximation (TFA). Numerical results, which 
allow us to outline existence regions of fundamental (single-peak) solitons combining the unstaggered and 
staggered components, and identify their stability (almost all the solitons are stable), are summarized in 
Section IV. Analytical results for the stability arc reported too, based on the Vakhitov-Kolokolov (VK) 
criterion for the two-component system. The numerical results corroborate the predictions of the VA quite 
well; in particular, it is confirmed that the mixed unstaggered-staggered solitons exist only in the case of the 
repulsive XPM, whilst the SPM is self-attractive in both components. The paper is concluded by Section V. 

II. THE MODEL 

The underlying system of the DNLS equations for lattice fields <p n and ip n is 

= ~ (4>n+l + 4>n-\ - ^4>n) ~ {\4>n\ 2 + ft |^n| 2 ) 0n> (!&) 

i 'dl^ n = _ 2~~~ + ^ n ~ 1 ~ 2 ^' n ) ~~ (l^«| 2 + P l^" ! 2 ) V'ra, ( lb ) 

where t is time in the case the BEC mixture, or the propagation distance in the array of optical waveguides, 
m is the relative atomic mass of the two species in the case of BEC, or the inverse ratio of the inter-site 
coupling constants in the waveguide array, and (3 is the relative coefficient of the XPM coupling between 
the fields, assuming that the coefficients of the self-attractive SPM nonlinearity for both fields are scaled to 
be 1. It should be mentioned that the model based on Eqs. ([1} is not the most general one, as, rescaling 
both fields to make their SPM coefficients equal to 1, one can make the XPM interaction asymmetric, with 
two different coefficients in Eqs. (pji) and (b), ^ fy. Nevertheless, quite generic results concerning the 
discrete solitons can be obtained within the framework of the present system. 

Solutions with unstaggered n and staggered ip n components and two chemical potentials, A and /i, are 
sought for as 



4> n {t) = e- lXt u n , Mt) = e^' (-1)% 



(2) 



where real u n and v n satisfy the following stationary equations, 



3 



(A - 1) U„ + - (u n +l + Un-i) + (u n + /3v n ) U n = 0, 
/" - — J V n - — (V n +1 + V n -l) + {v n + fan) v n = 0, 

to / 2m 



that can be derived from the Lagrangian, 
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In the large-m limit, which is tantamount to the TFA [Toj for discrete equation (j3]o) , this equation demon- 
strates that v n can be eliminated in favor of u n , hence in this case the coupled stationary system reduces to 
a single equation. 

In the next section, wc present variational solutions based on an exponential ansatz for fundamental 
(single-peak) solitons, and continue the analysis in Section III by means of numerical methods. For given 
(3 and to, we determine regions in the (A, fj,) plane for which single-peak numerical solutions exist and are 
stable. It is also found that the related energy surfaces, i.e., norms of the two components as functions of 
A and /i, always decrease in A and either increase or decrease monotonically in /i, depending on the sign 
of 1 — 1 . In this way, the generalized Vakhitov - Kolokolov (VK) stability criterion for two-component 
solitary waves can be applied here [TH, [l7| . In related two-component continuous systems fTsM2lT | , modeled 
by coupled continual NLS equation, one can introduce a new parameter (the ratio of A and /z) and rescale 
the variables, to make the stationary states depending on one (rather than two) effective chemical potential 
(2lj . Moreover, a generalized VK stability criteria was developed for a system of N incoherently coupled 
continuous NLS equations in Ref. [22| . 

As for discrete systems, the single DNLS equation with the arbitrary power-law nonlinearity was studied, 
by means of the variational approximation (VA), in Ref. (23| . and the stability of multi-soliton bound states 
in the DNLS equation with the cubic self- focusing nonlinearity was investigated in Ref. (24[. A complex 
version of the VA made it later possible to make predictions about collisions between moving lattice solitons 
in the same basic model [25jj. Another variational ansatz, relevant for DNLS solitons located on or anywhere 
between lattice cites, was elaborated in Re f. [2611 . The VA was further generalized for the DNLS equation 
with the cubic-quintic on-site nonlinearity [27|]. Very recently, the accuracy of the VA-based description of 
static discrete solitons and their stability, based on ansatze with different numbers of free parameters, was 
investigated in a rigorous form in Ref. [28[. As concerns discrete two-component systems, the VA was used 
for studying the spontaneous symmetry breaking in parallel DNLS lattices, linearly coupled at all sites [29| . 
or at a single site [3(J ■ 



III. ANALYTICAL APPROXIMATIONS 

A. The variational approximation for the discrete solitons 

To apply the VA to the solution of Eqs. (EL we employ the exponential ansatz that was earlier used in 
the framework of other models [Hj], [26|, [ffi[30j: 



(5) 



We find the decay rates of the wave forms in Eq. ([5]), p and q, not from the variational principle, but b 
requiring the ansatz to satisfy the linearized limit of Eqs. ([3]) at n — > ±00: 



P = In (l - A + v/-A(2-A)) , 

q = In (mfi — 1 + \J mfi (mfi — 2)^ . 

For p and q to be real and positive, the allowed ranges of chemical potentials [i and A are 

2 

A<0, \i >0. 
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Substituting ansatz (JS|) into Lagrangian and carrying out the summation yields the effective Lagrangian, 

2L cS = -A 2 tanh (-) + — tanh (-) + XA 2 cothp + ( fi - B 2 coth q 
\2J m \2J \ m J 



A A B 4 
+ — coth (2p) + — coth (2q) + f3A 2 B 2 coth (p + q), 

which gives rise to the variational equations, dL c g/d (A 2 ) = dL e g/d (B 2 ) = 0, i.e. 

A 2 coth (2p) + \p coth (p + q)] B 2 = tanh (|) - A cothp, 



[/3 coth (p + q)\ A 2 + B 2 coth (2q) 



-— tanh f-) — ( /i I coth g 
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As seen from Eq. (JSJd) and ([7]), solutions to the variational equations with positive A 2 and B 2 do not exist 
in the case of (3 > 0, but a solution may exist at j3 < 0. 

The fact that the fundamental solitons of the mixed unstaggercd-staggered type may exist as the bound 
state of two components, which, in isolation, support solely ordinary unstaggered solitons (throu gh t he self- 
attractive SPM), suggests to identify the solitons of the mixed type as symbiotic ones, cf. Ref. [3l|, where 
symbiotic solitons were defined in the opposite case, for the continual system with the self-repulsive SPM and 
attractive XPM nonlinearities. On the other hand, the staggering effectively reverses the signs of the SPM 
nonlincarity and external potential, therefore, in the presence of a large-amplitude unstaggered component, 
the staggered one may be considered as a soliton with the intrinsic self-repulsive nonlincarity, trapped in 
the attractive external potential. Such a mode tends to exist and be stable, unless the effective intrinsic 
self- repulsion is too strong, making the existence of the trapped mode impossible [32j . 

We also note (this remark will be relevant for comparison with some numerical results presented in the 
next section) that a solution to Eqs. ©, considered as a linear system for A 2 and B 2 , may not exist when 
the determinant of the system vanishes, i.e., 



coth(2p) coth(2o) - /3 2 coth 2 (p + q) = 0. 



Nevertheless, a solution is possible under condition (|TU| if the right-hand sides of Eqs. 
same way as the two rows of the degenerate determinant, i.e., 



(10) 



are related in the 
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B. Three-layer solitons for 1/m —¥ (the discrete Thomas- Fermi approximation) 
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There is another case in which we can determine properties of the solution in an analytical form. When 
the staggered species is very heavy, i.e., m — > oo in Eq. |T|d), the second equation from system ([3]), at lowest 
order, takes the local form: 

[fjL+(vl + Pu 2 n )] v n = 0. (12) 

Equation (|12[) has three possible solutions, viz. 

vl = -M - Pu 2 n , (13) 

or v n = 0, which may be used to eliminate v n in favor of u„, cf. a similar approach allowing one to eliminate 
a heavy fermionic component in Bose- Fermi mixture [32j. Accordingly, discrete solitons, composed of three 
layers, can be built as follows: in the central region (inner layer), we use relation (|13p and substitute it into 
the first equation of system ©, which yields 

[(A - 1) - fy] u n + ^ [u n+1 + u„_x) + (1 - 2 ) ul = 0, (14) 

i.e., the stationary DNLS equation which gives rise to soliton solutions. Requiring this solution, in the central 
region, to be a part of a discrete soliton with a single peak and centered at n = 0, then one must have /3 2 < 1 
and A < /3/j,. 

It follows from Eq. that, since \x > [see Eq. (J7J], for to be positive, one must take — 1 < f3 < 0, 
which yields v 2 . = (-/?) (u 2 n - U 2 ) , with U 2 = ~n/f3 > 0. Provided that > U 2 , one thus has i>l > at 
n = and in some region around n — (in the inner layer of the solution, as defined above), However, the 
positiveness of the so defined will be lost at |n| > with N large enough, as u 2 n for soliton solutions 
decays at n — > oo. Thus, for n > N and n < — N (in the two outer layers), the discrete mode can be 
extended upon taking the other root of Eq. (|12p for v n , namely, v„ = 0, which thus causes u„ to satisfy the 
usual DNLS equation, following from Eq. ([3]b) with v n = 0: 

1 3 
(A - 1) u n + - (u n+1 + u n -i) + u n = 0. (15) 

Obviously, Eq. (|T5|) has usual solution vanishing as \n\ —> oo for A < [recall A < is imposed by Eq. ([7])], 
thus the composite soliton can be constructed by combining the appropriate solutions in the inner and outer 
layers. The conditions of matching the discrete fields at n = ±N, which includes setting = (as required 
by the TEA), imposes two constraints on the set of parameters A, /i, and N, hence the solution is expected to 
exist along a curve in the plane of (A,/Lt), which is corroborated by numerical findings presented in the next 
section. Note that, in the framework of the present approximation, there is actually no difference between 
the unstaggered and staggered forms of the solution for v n , as only is determined by Eq. (fT5T) . 



IV. NUMERICAL SOLITON SOLUTIONS 

A. The formulation of the numerical problem 

We look for numerical solutions to Eqs. for spatially symmetric solitons, with it_ n = u +ni v-„ = 
and both fields u„ and v„ monotonously decaying with the increase of n, but never changing their signs, 
to support the unstaggered and staggered shapes of the underlying components 4> n and ip n , respectively, 



according to Eq. ©. At n = 0, Eqs. (0 yield 
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«i = - [u§ + pvl + (A - 1)] m , (16a) 



P u o + ( M — — 



wo- (16b) 



According to the above conditions, solutions to Eqs. (|16|) must satisfy constraints uo > «i > and i>o > 
«i > 0, thereby implying that 

- 1 < u 2 + pv 2 Q + (A - 1) < 0, (17a) 
1 2 

— < vl + fiul + fj, < — . (17b) 
m m 

Continuing in this manner, i.e., imposing bounds Ui > U2 > 0, v\ > V2 > 0, and so on, as it follows from 
Eqs. ^ at n = 1, 2, one successively restricts the region of the (A, /z) plane in which the soliton solutions 
are possible. 

The numerical solution of Eqs. ((3|) was carried out by means of a discrete version of the shooting 
method, which used the VA-prcdictcd solution as the initial guess, and was iterated until discrete wave forms 
monotonously decreasing with n without the change of the sign, up to the level of (u„, v n ) ~ 10 -5 (uq, vq), 
were found. 

For stability testing, we introduced initial perturbations, multiplying the stationary solutions by 

[1 + 5 exp(ifcn)], (18) 

with perturbation amplitude S ~ 5% and k N ~ 1, where N is the effective size of the discrete soliton. Then, 
the evolution of the thusly perturbed solution was simulated forward in time until t = 50. The results of the 
simulations were characterized by "stability numbers" S$ and for the two components, which are defined 
as root-mean-square changes in the relative amplitude of the solution, compared to the initial values, over 
the part of the lattice where the discrete soliton is located. For stable solutions, we obtain \S$,S^\ <C 1, 
while for unstable ones |S^, S^\ grow to values > 1. 

Since these solitons are symbiotic, one might suspect that they could be unstable to efforts to pull their two 
components apart. We have also checked for this possibility numerically, as above, by taking wavenumbers 
k u and k v with opposite signs in perturbation factors (fT5|) for the two fields. All solutions that we tested 
in this way, which had tested out to be stable against other perturbations, were found to be stable in this 
sense too. 



B. Dependence of solutions on the parameters 

In agreement with the prediction of the VA, numerical solutions for the solitons were found solely for 
j3 < 0, and, as suggested by Eq. (|14[) . /3 = — 1 is a critical value. When m = 1, we can demonstrate this 
with numerical results which makes it possible to identify two distinct cases, /3 < — 1 and — 1 < (3 < 0, seen 
in Fig. 1. When j3 approaches —1 from either side, we find solutions in a region which shrinks toward the 
line 

M = 2-A (19) 

in the (A, ^)-plane. It is worthy to note that, as can be found from inspection of Eqs. © and ([TTj) . both these 
equations reduce precisely to Eq. (fT9|) in the case of m = 1 and j3 = — 1, i.e., only the "double-degenerate" 
solution selected by Eqs. © and (fTTj) survives in this case. Note also that Eqs. © with m = 1 yield 
equal decay rates p and q for the two components of the soliton exactly under condition Eq. (fT9|). i.e., the 
soliton surviving in the limit of m = 1 and f3 = — 1 is characterized by equal localization lengths of the two 
components. 

Moving away from the critical value, j3 = —1, in either direction, Fig. 1 shows that the existence region 




(a) (b) (c) 

FIG. 1: (Color online) Existence regions (black) for the numerically found single-peak discrete solitons in the (A, /i) 
plane, for m = 1 and (a) /3 = —0.5, (b) ft — —1.01, (c) — —2. The dashed blue line corresponds to /i = 2 — A [see 
Eq. (|19l) ]. to which the existence region shrinks in the case of p = — 1. 

of the numerically found solitons widens, and, simultaneously, the region moves away from line (|19[) , staying 
on one side of this line, depending on the sign of 1 + /?. For /3 sufficiently far from —1, additional solution 
regions begin to split off from the primary one in the (A, fi) plane. Additional solution branches break off 
from the primary "trunk" at small |A|, then shrinking and disappearing as |A| increases, while the main 
trunk widens as A — > — oo. 

With the increase of the relative-mass parameter, to, the existence region of the soliton solutions in the 
(A, n) plane shrinks, following Eqs. (fT7| . This trend is observed in Fig. 2, which suggests that the region 
contracts toward a line in the (A, /x) plane at to — s- oo (as predicted by the TFA presented above). In case /3 is 
far enough from —1 to permit additional branches in the existence diagrams, we observe that such branches 
collapse into the primary one (the "trunk"), which then itself collapses into a line, as can be seen in Fig. 
2 for j3 = —2. It is also worthy to note that the bottom boundary of the existence region in Fig. 2 moves 
upward with the increase of to at fixed A. 

Examples of the solitons, including the juxtaposition of their numerically found and VA-predicted profiles, 
are displayed in Fig. 3 for /3 = —0.5 and m = 1. Fixing A = —4, we pick solutions from the larger lower 
stability region and the upper one in Fig. 1(a) corresponding to fi = 2.6 and (i = 3.55, respectively. Note 
that, while the profiles of the unstaggered component u n are very similar to one another, the solutions for 
v n are different. In both cases, the variational solutions agree well with their numerical counterparts. 

In Fig. 4, the solitons are plotted for the four branches of the existence region in Fig. 1(c) when j3 = — 2 
and m = 1. The first solution, shown in Fig. 4(a), belongs to a very narrow existence branch, which is barely 
discernible in Fig. 1(c) (its vertical width is A/i < 0.002), and exists along the bottom right of the main 
existence region (near the edge where A ~ — 1.3 and fi ~ 5.0). The other solutions are taken from the large 
lower existence region [Fig. 4(b)], the large upper one [Fig. 4(c)], and the thin upper stripe which splits off 
from the large upper branch [Fig. 4(d)]. Solutions from the lowest region [Fig. 4(a)] feature wider profiles 
in u n (note that both u±\ for them are on the same order of magnitude at uq) than do the solutions from 
all the other branches, which exhibit sharp profiles and agree well with the VA. On the contrary, the broad 
profile for u n in Fig. 4(a) cannot be approximated properly by the exponential ansatz ([5]). 

Because the numerical method employed here starts in a region where the variational equations, Eqs. 
have a solution, we cannot be absolutely sure that numerical solutions exist only in the dark areas shown 
in Figs. 1 and 2. In principle, other branches of numerical solutions might exist too, being unrelated to the 
VA, although this does not seem plausible. 

C. Soliton stability 

Systematic simulations of the evolution of perturbed solitons shown in Fig. 3 and Fig. 4 confirm that 
they are stable, with the exception of the one in Fig. 4(a). Further, systematic tests clearly suggest that the 





numerically found solitons are stable if their shapes are close to those predicted by the VA, whereas "broad" 
solutions, which disagree with the VA, turn out to be unstable. Actually, such unstable solitons are found 
only near the lower boundary of the regions shown in Fig. 1. 



n n 




(c) (d) 

FIG. 4: (Color online) The same as in Fig. 3, for f3 = -2, m = 1 and (A, fx) = (-1.25, 5.58) (a), (A, fi) = (-1.25, 6.005) 
(b), (A, fj.) = ( — 1.25, 9) (c), and (A, fi) = ( — 1.25, 10.7) (d). For the soliton in panel (a), the variational approximation 
provides a poor fit to the numerical solution, and this soliton is unstable. Other solitons are well approximated by 
the variational ansatz, and are stable. 

In order to deduce the stability in a more general way, we define the energies (norms) of the components, 

oo oo 

w u {\n)= K\\w v (\,n)= M 2 - ( 2 °) 

n— — oo n— — oo 

In Fig. 5 W U (X, fi) and W v (\, fi) are plotted for m = 1 and (3 = 0.5, and in Figs. 6 and 7 we do the same 
for f3 = —1.1 and /3 = — 2, respectively. There is a noticeable difference between the energy surfaces for the 
P < — 1 and — 1 < (3 < cases. In all cases considered, the energy surfaces are monotonous in A and fi 
over the stability regions (this finding agrees with the stability results reported in Ref. [23|). However, for 
— 1 < f3 < 0, the energy surfaces increase in A (holding /j, fixed) and decrease in fi (holding A fixed), as seen 
in Fig. 5. The opposite feature is observed at (3 < —1: the energy surfaces decrease in A at fixed fi, and 
increase in fi at fixed A. 

More can be stated about the stability by means of the VA. The substitution of ansatz ([S]) into Eqs. (|2T))) 
yields 

oo 

W u (X,fi) = A J2 e" 2pH = Acoth(p), (21) 




(a) (b) 

FIG. 5: (Color online) Color-coded plots of the energy surfaces W U (X, /i) (a) and W V (X, /i) (b) in the existence regions 
for the discrete solitons in the (A, fi) plane, at m = 1 and /3 = —0.5. The energy increases in A (holding fi fixed) and 
decreases in /j, (holding A fixed). 
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(a) (b) 

FIG. 6: (Color online) The same as in Fig. 5, but for m = 1 and /3 = —1. 01. The energy surface decreases in A 
(holding n fixed) and increase in /i (holding A fixed). 

oo 

W v {X,fJ,)=B e ~ 2<iH = Bcoth(q). (22) 

n— — oo 

where A and B are functions of A and \x determined by Eqs. (©). As is known from the generalized VK cri- 
terion for systems with two conserved norms [l?}, a stability change occurs when Jacobian d(W u , W v )/d(X, /i) 
changes its sign. We follow this approach in Fig. 8, where the zero locus of the Jacobian is plotted, along 
with the region of the existence of the numerically found solitons, for j3 = — 2 and m = 1. It is observed that 
the stability change predicted by the VA nearly coincides with the lower boundary of the existence region. 
The agreement is not perfect since the VA does not produce exact results, but the mismatch is quite small. 
The majority of the soliton solutions, which are located above the stability-change locus, are stable; unstable 
are the solitons, such as the broad one displayed in Fig. 4(a), which are found in a tiny area adjacent to the 
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FIG. 8: (Color online) The numerically found existence regions (black) in the (A, /x) plane for the single-peak solitons, 
juxtaposed with the (red online / dark grey in print) curves defined by d(W u , W v )/d(X, n) = 0, as produced by the 
variational approximation, for m = 1 and /? = —2. The solitons existing above the red line are stable. 

lower boundary which is actually bounded by the Jacobian's zero locus crossing the existence region. 

Comparing Figs. 5-7, we see that as j3 — > — 1, the energy values for the obtained solitons increase. We 
observe that this, in turn, corresponds to a change in the stability, and it was found that solitons occurring 
in the narrow region shown in Fig. 1(b) for f3 = —1.01 arc unstable. One such soliton is plotted in Fig. 
9. Notice that the unstable numerical solution is again essentially wider than its variational counterpart, 
i.e., as in Fig. 4(a), the variational approximation is a poor fit to the broad soliton. We have tested the 
stability of similar soliton solutions for (3 = —1.10 and fi = 6.0, and have found them to be stable. So the 
instability region appears to be localized around /3 = — 1. Lastly, direct simulations demonstrate that the 
unstable broad soliton solutions, such as the one displayed in Fig. 4(a), decay into a combination of multiple 
breathers and emitted radiation, as seen in the example in Fig. 10. 
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FIG. 9: (Color online) An example of a discrete soliton for /3 — —1.01, m = 1 and (A,/i) = (—3.55, 6). Symbols and 
lines depict the numerical solutions and prediction of the variational approximation, respectively. The soliton shown 
here is unstable. 
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FIG. 10: (Color online) Plots of (a) <j>„(t) and (b) ip n (t) of the evolution of the unstable soliton from Fig. 4(a), which 
corresponds to /3 = —2, m = 1, (A, /i) = ( — 1.25,5.58). Note the escaping radiation in both components, and the 
oscillating breathers which are left behind. 

V. CONCLUSIONS 

We have introduced the symmetric system of DNLS (discrete nonlinear Schrodingcr) equations with the 
self- attractive on-site SPM nonlinearity and repulsive XPM interaction, which supports two-component 
solitons of the symbiotic unstaggcrcd-staggcrcd type. The system may be implemented in a mixture of 
two BEC species with identical or different atomic masses, and (in principle) in arrays of bimodal optical 
waveguides. In the analytical part of the work, the VA (variational approximation) was developed, based on 
the exponential ansatz for the fundamental (single-peak) solitons. In the limit of the large relative mass of the 
two species, the TFA (Thomas-Fermi approximation) was elaborated too, which reduces the coupled system 
to two different single-component DNLS equations in the inner and outer layers of the solutions. Further, 
by means of the numerical solution we have identified areas in the plane of the two chemical potentials 
(propagation constants) where discrete solitons exist. It has been inferred that the VA and TFA agree well 
with the numerical solutions, except for a stripe near the lower existence boundary, where broad solitons are 
poorly approximated by the exponential ansatz. Direct simulations of the evolution of the perturbed solitons 
demonstrate that all the solitons which are well approximated by the VA (i.e., almost all the solutions) are 
stable. Only the broad solitons, which are not accommodated by the VA, are unstable. The results for the 
stability can be accurately predicted by means of the generalized VK criterion for the two-component system 
(with the stability change corresponding to the vanishing of the respective Jacobian), realized in terms of 
the VA. 



It may be interesting to extend the work by considering multi-soliton (multi-peak) bound states of tn§ 
unstaggered- staggered type. A challenging problem is to generalize the system for two-dimensional lattices 
and various types of discrete two-dimensional solitons, including solitary vortices. 
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